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Abstract 



Multi-scale wave propagation problems are computationally costly to solve 
by traditional techniques because the smallest scales must be represented over 
a domain determined by the largest scales of the problem. We have devel- 
oped and analyzed new numerical methods for multi-scale wave propagation 
in the framework of heterogeneous multi-scale method. The numerical meth- 
ods couples simulations on macro- and micro-scales for problems with rapidly 
oscillating coefficients. We show that the complexity of the new method is 
significantly lower than that of traditional techniques with a computational 
cost that is essentially independent of the micro-scale. A convergence proof is 
given and numerical results are presented for periodic problems in one, two 
and three dimensions. The method is also successfully applied to non-periodic 
problems and for long time integration where dispersive effects occur. 



1 Introduction 

We consider the initial boundary value problem for the scalar wave equation, 



on a smooth domain Q c R N with A e (x) a symmetric uniformly positive definite 
matrix. We assume thatA E has oscillations on a scale proportional to e <C 1. The so- 
lution of ([I]) will then also be highly oscillating in both time and spatial directions 
on the scale e. It is typically very computationally costly to solve these kinds of 
multi-scale problems by traditional numerical techniques. The smallest scale must 
be well represented over a domain, which is determined by the largest scales. For 
wave propagation small scales may also originate from high frequencies in initial 
data or boundary data. We will however focus on the case when they come from 
strong variations in the wave velocity field. Such variable velocity problems occur 
for example in seismic wave propagation in subsurface domains with inhomoge- 
neous material properties and microwave propagation in complex geometries. 



u s - 

tt 

u e = 



■V-A E Vu £ = 0, nx{0<t<T}, 

f, u\ = g, nx{t = o], 



(i) 



Recently, new frameworks for numerical multi-scale methods have been pro- 
posed, including the heterogeneous multi-scale method (HMM) JSJ and the equa- 
tion free methods II 1 311 - These methods couple simulations on macro- and micro- 
scales. We use HMM, [5|[6j|4l, in which a numerical macro-scale method gets 
necessary information from micro-scale models that are only solved on small sub 
domains. This framework has been applied to a number multi-scale problems, for 
example, ODEs with multiple time scales II10II . elliptic and parabolic equations 
with multi-scale coefficients [0 [T21 HI, kinetic schemes [6] and large scale MD 
simulation of gas dynamics II15II . 

On the macro-scale we will assume a simple flux from, 

u tt -V-F = 0, (2) 

in our HMM approximation of the wave equation ([I]) . The solution u should be a 
good approximation of the solution to ([!]) and the value of F on the macro-scale 
grid is computed by numerically approximating on small micro-scale domains. 

The goal of our research is to better understand the HMM process with wave 
propagation as example and also to derive computational techniques for future 
practical wave equation applications. One contribution is a convergence proof in 
the multidimensional case that includes a discussion on computational complex- 
ity. The analysis is partially based on the mathematical homogenization theory for 
coefficients A e with periodic oscillations ]2l [3] • 

Classical homogenization considers partial differential equations with rapidly 
oscillating coefficients. As the period of the coefficients in the PDE goes to zero, 
the solution approaches the solution to another PDE, a homogenized PDE. The 
coefficients in the homogenized PDE has no e dependency. For example, in the set- 
ting of composite materials consisting of two or more mixed constituents (i.e., thin 
laminated layers e periodic), homogenization theory gives the macroscopic prop- 
erties of the composite. It is an interesting remark that the macroscopic properties 
are often different than the average of the individual constituents that makes up 
the composite Q3J. The wave equation ([!]), withj4 E (x) = A(x,x/e) and A(x,y) is 
periodic in y, have an homogenized equation, 

a n - v -Avu = o, nx{o<t<T}, 

u=f, u t = g, Ox(t = 0), 

where A(x) is called the homogenized or effective coefficient. The homogenized 
solution u can be used as an approximation of the solution u s of the full equa- 
tion since i/(x) = u(x) + 0(e). Note that, the homogenized equations are often 
less expensive to solve with numerical methods, since the coefficients varies slowly 
without e variations. We refer to [2j[T8l[3l [THESIS ^ or more about homogeniza- 
tion in general. 

It should be noted that even if our numerical methods use ideas from homog- 
enization theory they do not solve the homogenized equations directly. The goal 
is to develop computational techniques that can be used when there is no known 
homogenized equation available. In the research presented here many of the ho- 
mogenized equations are actually available and could in practice be numerically 
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directly approximated. We have chosen this case in order to be able to develop a 
rigorous convergence analysis and to have a well-understood environment for nu- 
merical tests. We also apply the techniques to problems that does not fit the theory. 
In example 4.2.3 an equation with non-periodic coefficients is approximated and 
in example [4. 5 1 an equation is solved over very long time. The latter is particularly 
interesting since the homogenized solution contains dispersive effects, which in- 
fluence the solution for t = 0{e~ 2 ). This dispersive process is captured by a high 
accuracy HMM technique without explicit approximation of any dispersive term. 

The article is organized as follows: In section[2]we discuss first the HMM frame- 
work in a general setting and thereafter in section |2.1| our HMM method for the 
wave equation. We give a rigorous proof of the approximation error by the HMM 
method in the periodic coefficient case in section[3] In section|4]we show numerical 
results, which also includes a non-periodic problem and an example with very long 
time. The last section [5] ends this paper with our conclusions. 



2 Heterogeneous multi-scale methods (HMM) 

In the HMM framework, the general setting of a multi-scale problem is the follow- 
ing: We assume that there exists two models, a micro model /(u, d) = describing 
the full problem and a coarse macro model F[u, d) = 0. The micro model is accu- 
rate but is expensive to compute by traditional methods. The macro model give a 
coarse scale or low frequency solution u, assumed to be a good approximation of 
the micro-scale solution u and is less expensive to compute. The model is however 
incomplete in some sense and requires additional data. We assume that F(u, d) = 
can still be discretized by a numerical method, called the macro solver. A key idea 
in the HMM method is to provide the missing data in the macro model (d) using a 
local solution to the micro model. The micro model solution iz is computed locally 
on a small domain with size proportional to the micro-scale. The initial data and 
boundary conditions (d) for this computation is constrained by the macro-scale 
solution u. 



2.1 HMM for the wave equation 



We will formulate a general HMM framework for the wave equation on the domain 
Y = [0, l] d . Let u e be Y-periodic and solving, 



u e n = V -A e Vu e , 7x{0<t<T(, 
u e =f, u e = g, yxft = 0|. 



(4) 



We follow the same strategy as in [1] for parabolic equations and in II19II for the 
one-dimensional advection equation. See also JS] . We assume there exists a macro- 
scale PDE of the form 



u tt - V-F(x,u,Vu,...) = 0, Y x {0< t < T}, 
"=/. u t = g, Yx{t = 0}, 

u, Y -periodic. 



(5) 
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where F is a function of x, u and higher derivatives of u. The assumption on ([5]) 
is that iz as u e when e is small. In the clean homogenization case we would have 
F = AVu, but we will not assume knowledge of a homogenized equation. Instead 
we will solve the PDE Q, only in a small time and space box, and from that 
solution extract a value for F. The form of the initial data for this micro problem 
will be determined from the local behavior of u. In the method we suppose that 
F = F(x,Vu). 

Step 1: Macro model discretization. We discretize ([5]) using central differences 
with time step K and spatial grid size H in all directions, 

' U n+l = 2u n_ u n-l + Kil F W _ p (X) )+...+ £ f F C<Q _ pW 



H \ m+-e! m--e 1 J H \ m+-e d m—-e d 

F n , = F(x m i, ,P n . ), k=l,...,d, (Note: F n , is a vector.) 

m-ie t rn--e k > m -l e /> m-±e t 

(6) 

where F^ ±1 / 2 e k * s ^ evaluated at point x m±1 / 2ei - The quantity P^ ± ^ ^ approximates 
Viz in the point x m±1/2 - We show an example in Figure|2.1|of the numerical scheme 



in two dimensions. There P^ + i e is given by the expression ( [7T| ) in the Appendix. 
Step 2: Micro problem. The evaluation of F n , in each grid point is done by 

m--e k 

solving a micro problem to fill in the missing data in the macro model. Given the 
parameters x m _i e and P n l , we solve a corresponding micro problem over a 

small micro box Y e , centered around x m _i e . In order to simplify the notation, we 
make a change of variables x — x i — > x. This implies that A e (x) — > A E (x + 
x m _i et ). The micro problem has the form, 

r u E n -V -A e Vu e = 0, Y £ x{0<t<T}, 

u E (x,0) = (P" , )-x, iz E (x,0) = 0, 7 £ x{t = 0}, (7) 

m--e k 

u e — iz e (x, 0), y £ -periodic. 

We keep the micro box size of order e, i.e. t, diam7 E = 0{e). We note that the 
solution u e is an even function with respect to t (i.e. u e (x, — t) = u E (x, t)) due to 
the initial condition u E (x,0) = 0. 

Step 3: Reconstruction step. After we have solved for for u e for all Y e x [0, t] 
we approximate F n , fa F(x m _i„ ,P" . ). The function F is the mean value of 

m--e k m z e k m--e k 

f E = A e Vu e over [— rj, rj] d x [— t,t] where [— 77, 77] d c Y e . The approximation 
can be improved with respect to the size of z/e and 17/e, by computing a weighted 
average of f e . We consider kernels K described in 1 1 1 Oil : We let K p ' q denote the 
kernel space of functions K such that K e C?(R) with supp K = [—1,1] and 



c 

K(t)t r dt 



1, r = 0; 
0, l<r<p. 
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Figure 1: The Numerical scheme ( [71] ) for P in two dimensions. In the figure above 
the two components of F in two different positions are given by F i+l / 2 j and G i J+1 / 2 - 
The U points involved in computing F n , = G t j+1/2 an d vu w P" , are indi- 

cated by filled circles. 
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Furthermore we will denote K„ as a scaling of K 

K n (x):=-K[-), K e K p ' 9 , 
17 yr]J 

with compact support in [—17,17]. We use kernels of this sort to improve the ap- 
proximation quality for the mean value computation, 



K,(t)K n (x)f°dxdt, f e =A e {x + x m+1 _ ek )Vu e , (8) 



where here the multi variable kernel K„(x) is defined as 

K 11 (x) = K n (x 1 )K v (x 2 )---K v (x d ), (9) 

using the single valued kernel K-, still denoted by K~. The domain Y e is chosen 
such that [— 17,17 ] d c Y e and sufficiently large for information not to propagate 
into the region [— 17,17] d . Typically we use 

y = [-y max ,y ma J d , y m ax = i7 + T\/sup|[A E || 2 , (10) 

c.f. discussion about micro solver boundary conditions in H19II . In this way we do 
not need to worry about the effects of boundary conditions. Note therefore that 
other types of boundary conditions could also be used in ([7]). 



Remark. It is possible to find functions with infinite q. In M0\l a kernel K exp is given, 
where p = 1 and q is infinite: 



[0, M > 1, 



where C is chosen such that f K exp (x)dx = 1. This kernel is suitable for problems 
where A 6 is of the form A e (x) =A{x/e). 

Remark. The weighted integrals above are computed numerically with a simple trape- 
zoidal rule. 

Remark. In our implementation, the micro problem ([7]) is solved with the same nu- 
merical scheme as the macro problem ([6]). 

2.2 Computational cost 

Let us assume that the time step is proportional to e in all direct solvers. Using a 
direct solver for ((4]) on the full domain implies a cost of order e~ (d+1 \ The total 
cost for HMM is of the form (cost of micro problem) x M d where M d is the number 
of micro problems needed to be solved per macro time step. The cost of a single 
micro problem is of the form (t/e) x (tj/c) 1 *. We assume kernels with t,?? ~ e 
and that M d does not depend on e. With these assumption our HMM method 



6 



has a computational cost independent of e. The constant can, however, still be 
large. Fortunately the computational cost of the HMM process can be reduced 
significantly. We observe that the function §8§ is linear in p. It is in fact composed 
of three linear operations: 

1. Compute initial data u(x, 0) and u t (x, 0) from p, u(x, 0) = p ■ x. 

2. Solve u e n - V • A £ Vu e = for < t < t. 

3. Compute average F = JJ K T K n f e dx dt where f e =A e Vu e 

The first operation is clearly a linear operation. In step two we compute a solution 
to a linear PDE, therefore this step is linear as well. Computing the integral average 
in step three is also a linear operation. 

As a corollary we can apply the HMM process to a smaller number of micro 
problems and form linear combinations of those for any given F computation. More 
precisely, after precomputing F(x,e t ), i = 1,2,..., d we can compute F for fixed 
x e Q and anyp e R d , 

d 

F(x,p) = Y i PiF(x,e i ), (11) 

i=l 

where p ( is the ith coefficient in p in the basis e 1 ,e 2 , ■ ■ ■ ,e d . In conclusion, by 
precomputing the micro problems F(x m ,e;) in ( [TT] ) we only need to solve d micro 
problems in each macro grid point x m = mH. There is no need to solve any micro 
problems again in the next macro time step. The complexity is as before 0{X), but 
with a lower constant not depending on the number of time step. 

Remark. In fact, ifA e is e-periodic and the macro grid is such that x m = r( mod e\ 
where r is constant and independent of m, we only need to solve d micro problems in 
total. In this case, the total cost is independent of both e and the macro grid size H. 



3 Convergence theory 

In this section we apply the HMM process to the problem with A e (x) =A[x/e) 
where A is a Y -periodic symmetric positive matrix and show that it generates results 
close to a direct discretization of the homogenized equation ([3]). In particular we 
show that 

F(x,p) = F(i,p) + ^^j J. (12) 

The function F and F are defined in ([8]) and ([5]) respectively and we note that 
here F(x,p) =Ap. The integer q depends on the smoothness of the kernel used to 
compute the weighted average of f s in ([§]). 

We will formulate the problem in the setting of elliptic operators. For the anal- 
ysis we solve the micro problem ([7]) over all of R d 

u e t - V-A e Vu E = 0, R d x{0<t<T}, 
u > (13) 

u e =p-x, u e = 0, K d x{t = 0}. 
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Note that this gives the same F as in ([8]) if we choose a sufficiently large box Y e . 

Theorem 1. Let F(x ,p) be defined by §8§ where u e solves the micro problem ( |13| ), 
A e (x) = A(x/e) and A is Y-periodic and smooth. Moreover suppose K e K p ,q , f and 
g is smooth and t = r\. Then for p / 0, 

^\P(x ,p)-F(x ,p)\<c(^\ , 

where C is independent of e, r\, p and q. Furthermore, for the numerical approxima- 
tion given in ([6]) in one dimension, with H = ne for some integer n and smooth initial 
data, we have the error estimate 

\Ul - u(x m , t„)| < C(.T) [H 2 + (s/tj)«) , < t B < T, 

where u is the homogenized solution to 

Proof. We will prove the Theorem in the following steps: 

1. Reformulate the problem as a PDE for a periodic function. 

2. Define an elliptic operator L(y). 

3. Expand V y • A(y) and v(y, t) (to be defined) in eigenfunctions to L(y). 

4. Compute time dependent v,(t) coefficients in the above eigenfunction expan- 
sion. 

5. Compute the integral of f e to get F. 

6. Compute the solution to a cell problem and give final estimate. 

Step 1: Express the solution to ( [13] ) as 

u £ (t,x)=p-x + v(x/e, t). (14) 

We insert this into ( [131 ) to get a PDE for v 

v„ = \V y -A{y)p + iV y ■A(y)V y v(y), 
v(x,0) = 0, v t (x,0) = 5 

where y = x/e. Since A is Y -periodic, so is v, and we can solve ( [151 ) as a Y -periodic 
problem. 

Step 2: We define the linear operator L(y) := — V y • A(y)V y on Y with periodic 
boundary conditions. Denote by Wj[y) the eigenfunctions and Pij the correspond- 
ing (non-negative) eigenvalue of L. Since L is uniformly elliptic, standard theory 
on periodic elliptic operators informs us that all eigenvalues are strictly positive, 
bounded away from zero, except for the single zero eigenvalue ||14ll 

= A < Aj < X z < ■■■ (16) 

and Wj e C°° forms an orthonormal basis for L 2 er (Y). Note also that w = l^l -1 is 
a constant function. 
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Step 3: We express V y • A(y) and v(y) in eigenfunctions of L: 

00 CO 

V / A(j) = ^a J w J (y) and v(y, t) = ^v/Ow/y). (17) 

Note that here a ; are column vectors and as in the one dimensional case we have 
that a Q = since the mean value of V y -A(y) is zero, 

V J -A(y)dy = 0. (18) 



y 
c 



V y -A(y)w (y)dy 

F 



1^1 



Step 4: We plug the eigenfunction expansions (17) into fll5| ) and find that 



Y i v^ , w j =-p-Y i a j w j --^Y i Lv j w j = J]—^w j -Y l ^ j w j . (19) 

;=0 j=l ;=1 ; = 1 j=l 

By collecting terms of Wj we get 

v' + 4v ; . = ^. (20) 

This is a system of ODE:s similar to the form, 

y" + ay = /3, (21) 
which has the solution of the form (a > 0) 

y{t)=Ae it ^ + Be~ it ^+-. (22) 

a 

Note that all Xi > ( j > 0) so it is known that the Vj functions in the problem have 
the form, 

v(t) = Ae — + B.e + r., r. = — (23) 

and the special v is given by 

P • Q n i 

v (t)= r + Ct + D = Ct + D since a = 0. (24) 

2e 

By plugging the general solution ( [23] ) into the initial conditions of ([15]), we can 
formulate equations forAj and Bj (j > 0), 

00 

v(0, x) = => V Vj(0)w%x) = => Vj(0) = => Aj + B } + r } = 0; (25) 

v t (0, x) = => J] v;(0)wj(x) = => v'(0) = => -?-J. Aj - JLJ -B } = 0. (26) 
j=o e £ 



Similary, for v (t) 

v (0) = =^> C = 0, and v' Q (Q) = =^> D = 0, 
thus v (t) = 0. We solve for A: and B ; and get 



r ■ ep • a ■ 
A i =B, = - = — -, =1,2,... 



All in all, the v ; (t) coefficients in explicit form are 



v (t) =0, 



2A, 



e « + e > 



?0 



£p-a, ep-a, I tJX, 

+ = I 1 - cos 



The solution to our problem ( [151 ) can then be expressed as 



V 



cos 



Step 5: Now plug the expression fl30"] ) into the expression ( fT4l ) 



(27) 



(28) 



J'=l,2,... 
(29) 

(30) 



V y -w/x/e) I A(x/e). 

(31) 



We write down and analyze the function f e in two parts f e =p- (A x + A 2 ), where 



A 2 (x/e, t) = -1^ I cos ^V y • Wj (x/eWx/e). 



(32) 



Step 6a: First we show that A x = A. To do that we need to use the so-called cell 
problem (or corrector problem, see Section 4.5.4 in llllll ). 



(L(y) X = -V y -A, Y, 

[X Y -periodic. 

We rewrite the cell problem ( [33) ) using a eigenfunctions expansion 

OO 00 _ 

2 ^jXjWj = ~Yi a i w j ^Xj = -y, J = 1,2, ... , 

;=0 j=l A J 



(33) 



(34) 
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where Xj are column vectors with coefficients of % in the eigenfunctions expansion. 
For the other term A 1 we now will make good use of the eigenfunction expansion 
of the cell solution %, 



K,(t)K v (x)A 1 dxdt= | K„(x) ^ + ^^-V y w J (x/e)jA(x/e)dx 



K v (x) (/ - V yX (x/e))A(x/e)dx 
K^x) (A(x/e) - V yX (x/e)A(x/s)) dx 
K n {x)(Atxle)-Kxle)V yX {xlej) dx 
A+ 6 



(35) 



where we used Lemma [T] in each coordinate direction. 
Step 6b: Now we should show that 



K^t)K (x-x )A 2 dt — 0, -->oo. 
' s 



For that we need a Lemma from [flO 



(36) 



Lemma 1. Let / e (t) = /(t, t/e\ where /(t,s) is 1-periodic in the second variable 
and d r f{t,s}/dt r is continuous for r = 0, 1, ...,p — 1. For any K e KP' q there exists 
constants C 1 and C2, independent of e and r\, such that 



E = \K 71 *f £ [t)-f(t)\<C 1 r 1 P + C 2 



r 1 



fit): 



f(t,s)ds. 



If f = /(t/e) then we can take C x = 0. Furthermore, the error is minimized if 17 is 
chosen to scale with e q ^ p+q \ 

We now apply Lemma [l] to obtain 



K^t)cos^—dt 



< Co 



2ne 



■ C- 



q/2 



£^9 



(37) 



Let bj and the column vector g(y) be defined as 



if T (t)cos dt, g{y) = ^bjXjWjty), 



(38) 
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where we again used the solution to the cell problem ( |33] ) in the formulation of 
g(y). We then express J J K x K r .A 2 dx dt using g, followed by a change of variables: 



K T (t)K n (x - x )A 2 dx dt = - 



K v V y ■ g(x/e)A(x/e)dx 



r\x + 



17 



-K{x)A 



e 

r]x + x 



V-g 



B 

tjx + x 



dx 



dx. 



By doing integration by parts, using K[e k ) = K(—e k ) = (fc = 1,2, . . . , d), together 
with Cauchy-Schwartz inequality, we obtain 



-K(x)A 



tjx + x 



t?x + x c 



dx 



-i,i] d 



< 



\l-iM d 



0(1) column vector 



V • -K(x)A 



7]X + X 



dx 



2 fr)x + x 



dx 



-i,i] d 



1/2 



(39) 



which is bounded by C||g||( £ 2 « where C is independent of e, and 17. Finally we 
need to show that ||g|| — » 0. This is done by observing that 



J=l 



max Ha "(l 2 Y 



where |b max | is bounded by, 



— f-Y 



(40) 



(41) 



following the computations in ([37]). Then finally, we add our results from the 
calculations above and get, 



F[x ,p) = p- 



K z (t)K^x-x )f e dxdt 



■P- 
p ■ \ A+ 6 



K z (t)K n (x - x )(A 1 (t) + A 2 (x/e, t)) dx dt (42) 



This proves the Theorem. 
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Final step: Now we show the error estimate \U^-u{x m ,t n )\ < C(T)(H 2 +(e/T]) 9 ). 
We observe that F in the Theorem is of the form 



where A is e-periodic. By 



F = A(x)p 



|A(x)-A|<C - 
.V 



(43) 



(44) 



By choosing H = ne for some integer n, we find that the macro scheme ([6]) is a 
standard second order discretization of the problem 



u tt -A\0)u xx = 0, nx{0<t<T}, 
u(0,x)=/(x), u t = g, nx{t = 0}, 



(45) 



since A(x m ) = A(mne) = A(0) for all m. Hence, if g = (the result is true also for 

U m=2 ( /(Xm ~ y/M®tn) +fi*m + V^Wt„)) + 0{H 2 ). (46) 

On the other hand, the solution of the homogenized ((3]) with g = is 

u(.x m , t n ) = - (f(.x m - \fAt n )+f(.x m + \fAtS) ■ 
Therefore we get the error estimate 

/(x+ v / A(0)t)-/(x+ VAt) +C[T)H 2 



(47) 



\U^-u[x m ,t n )\< sup 

|t|<T 



<|/'Lr|VA- Va| + c(t)h 2 

<C(r)(H 2 + ( e /n)"), 
for < t n < T. This proves the Theorem. 



(48) 

(49) 
(50) 

□ 



4 Numerical results 

In this section we show numerical results when applying the HMM process to var- 
ious problems in one, two and three dimensions. The notation in the experiments 
in the d-dimensional setting (d = 1,2,3) is the following: We let Y = [0, l] d de- 
note the macro domain and e be the micro problem scale. We denote by H and K 
the macro grid size and time step respectively and for the micro-scale we denote 
by h and k the grid size and time step respectively We use explicit second order 
accurate finite difference schemes (see Appendix). 
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Figure 2: Convergence results, error \F — F\ plotted against -q/e (t = 17) for fixed 
e = 0.01 and where A e = A 1 with only fast scales. The dashed line corresponds to 
the (e/rj) 9 term in Theorem[T] The bottom figure shows results for the exponential 
kernel (see Remark[2TT]) and indicates super algebraic convergence rate. 



4.1 Convergence study of different kernels 

In Figure [2] and Figure [3] we present convergence results for the flux F in terms of 
rj/e. We use different type of kernels for the problem ( [51] ) withj4 £ (x) = A 1 (x/e) 
and A e (x) = A 2 {x,x/e) where A^y) = 1.1 + sin(27iy) and A 2 {x,y) = 1.1 + 
i(sin2n;x + sin27iy). We compare our numerical results to the theoretical bounds 
in Theorem [T] On problems with both fast and slow scales which is not directly 
covered by Theorem[l] we see a (slow) growth of the error as t, 17 — » 00 consistent 
with the general approximation result in Lemma[T] We plot (e/17) 9 and tj p separate 
with dashed lines. 



4.2 ID results 

The general form for the one-dimensional examples is: 

' t = d x A e u £ x , Yx{0<t<T}, 



(51) 

Uj=0, 7x{t = 0} 5 
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10° 



p=1, q=lnf 



10 



■10 



10° 



10' 



rj/e 



Figure 3: Convergence results, error \F — F\ plotted against r\je (t = tj) for fixed 
e = 0.01 and where A e = A 2 with both fast and slow scales. The dashed line with 
negative slope corresponds to the theoretical bound from the first term in Lemma 
1 and the dashed line with positive slope corresponds to the if term. 
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where Y = [0, 1]. We show some dynamics in Figure [4] where we solved ( [51] ) for 
the A 8 and / given in example one below. The homogenized solution to ( |5"T] ) will 
be of the form 

{ut^djai,, 7xf0<t<r}, 
u = f, u t = 0, Yx{t = 0}, 

where A is given by the harmonic average of A(x,y) over one Y-period, 

r 1 



o A(x,y)' 



(53) 



and x being held fixed. 



4.2.1 Example one 

The first wave propagation problem we choose A e and / as 

A 8 (x) = A(x/e), A(y)= 1.1 + sin27rj, 
/(x) = exp(-(x-x ) 2 /cr 2 ), x = 0.5, cj = 0.1. 

We can compute A from ([53]) with techniques from complex analysis 



(54) 



, 21 

A= \ = 0.458257569495584... (55) 

u 100 

We will solve ( [52] ) with a fully resolved discretization or direct numerical sim- 
ulation (DNS), discretized homogenized solution (HOM) and our HMM method 
(HMM). We have used s = 0.01, 17 = lOe. In Figure[5]we show a snapshot of the 
solutions these methods after time T = 1. We use a kernel (same in both time 
and space) K e K 5 ' 6 , that is K has 5 zero moments and is 6 times continuously 
differentiable. 

4.2.2 Example two 

We now consider a variation of (TsTI) where A 8 is defined as 



The homogenized operator A will not be constant but a function with explicit x 
dependence. We compute analytically A(x) to be 



A(x) = \J ct(x) 2 - p 2 a(x)= 1.1 + -cos27ix, /3 = -. 



(57) 



For this experiment we use e = 0.01, K = 2H, H = 3.33 • 10~ 3 . For the micro 
problem we use k/h = 0.5 and h = e/64. The kernel from IK 9 ' 9 . The small H is 
to lessen the effect of the numerical dispersion. We show results from T = 1 in 
Figure [6] 
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T=0. 000000 



T=0.295444 







0.5 





























0.5 1 
T=0.886332 














<•• ■ " 















1 




0.5 

T=1 .477220 




Figure 4: The dynamics of the problem ( |5T) ) for 4 snapshots at t, = 3/{i\/ A), < 
i < 3. Note the small oscillations superimposed on the smooth profile. We observe 
how the initial pulse separates in one left going and one right going pulse. The 
effect of the periodic boundary condition can seen as the waves pass each other at 
the boundaries between frames 2 and 3. 
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T=1 .000000 



T=1. 000000 




0.5 1 0.8 0.85 0.9 



Figure 5: A snapshot of two super imposed solutions to ( |5TT ) together with a 
zoomed section. 
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T=1 .000000 



1 

0.8 
0.6 
0.4 
0.2 

0; 





X 


DNS 
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0.5 
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0.5 1 

Figure 7: A snapshot of a direction solution to ( [56] ) and the HMM solution (left) 
together with the material coefficients from example 4.2.3 (right). 

4.2.3 Example three 

In the last one-dimensional example the macro equation is unknown, i.e. homoge- 
nization does not provide A. We define A e as a sum of many micro-scale oscillations 



A 



(x)=l.H — y sin27i — , e i 
^ p. 



90 + 5(i- 1)' 



(58) 



/(x) = exp(-(x-x )7c7 2 ), x = 0.5, a = 0.1. 



A plot of A s is shown in Figure [7] The numerical parameters for the macro-solver 
(HMM and homogenized) are H = 3.33 • 10" 3 , K = 0.5H. The micro solver uses 
t = 10e 3 , 17 = £ 3 , h = £ 3 /64 and k = 0.5h. The kernel K used, for both time and 
space, is K e IK 5,6 . The results are shown in Figure [7] 

4.3 2D results 

In this section we present the numerical results for a two dimensional wave propa- 
gation problem over the unit square Y = [0, 1] x [0,1]. 

4.3.1 Example four 

We define A e (x) by the diagonal matrix, 



A E (x) = diag(a E (x),a s (x)) 
a e (x) = a(x/e), a(y) = 1.1 + sin 271;^. 

The corresponding homogenized matrix A in ([3]), 



(59) 



A=diag(V0.21,l.l), 



(60) 



19 



1 




o 1 ' 1 1 1 1 

0.2 0.4 0.6 0.8 1 



Figure 8: Full numerical simulation when A e has only a fast scale. 

and as in ID the initial data / is defined as a Gaussian, 

/(x) = exp(-||x-x ||2/a 2 ), 
x =[0.5 0.5], cr = 0.1. 

We use the exponential kernel K exp e IK 1 ' 00 . We let T = 1 and the scale parameter 
e is set to 0.01. The macro scheme uses H = 3.33 • 10~ 3 and K = 0.5H. The micro 
scheme uses h = e/64 and fc = 0.5h. We show the numerical results in Figure [8] [9] 
and [H 



4.3.2 Example five 

We let A e (x) be defined by the diagonal matrix, 
U E (x) = diag(a E (x),a e (x)) 

[a e (x) = a(x,x/e), a(x,y) = 1.1 + i(sin27ix 1 + sin2n;y 1 ). 
and the corresponding homogenized matrix A in ([3]), 



(62) 



|A(x) = diag(d(x),l.l), 

la(x)= V 'a(x) 2 - /3 2 , a(x) = 1.1 + 0.5 sin27ix 1 , /3 = 0.5. 



(63) 



The numerical parameters are chosen the same as in example |4.3.1 We show the 
numerical results in Figures [TT] and 12 
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0.2 0.4 0.6 0.8 1 



Figure 9: Direct solution of the homogenized equation when A is constant. 




i , , , , 1 

0.2 0.4 0.6 0.8 1 



Figure 10: HMM approach, when A e has only fast scales. 
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4.4 3D results 



Here we present numerical results for a wave propagation problem in three dimen- 
sions in a locally periodic media over the box Y = [0, l] 3 . 

4.4.1 Example six 

In this three dimensional problem A e [x) is a diagonal matrix 



,¥(x) = diag(a E (x),a E (x),a £ (x)), 
a e (x) = aO/e), a(y) = 1.1 + sin2ny 1) 



/(x) = exp(-||x - x \\l/a 2 ), 
x =[0.5 0.5 0.5], ct = 0.1. 



(64) 



and the corresponding homogenized matrix A in ([3]) is 

A(x) = diag (V0.21, 1.1, l.l) , (65) 
and the initial data / is defined as a Gaussian, 



(66) 



In this experiment we have used e = 0.01. The homogenized simulation uses 
T = 0.25, H = 0.05, K = 0.25H. The HMM solver uses T = 0.25, H = 0.05, 
K = 0.25H on the macro solver. The micro solver uses 17 = e, t = 5e, h = e/ '64, 
fc = 0.3h and a polynomial kernel K e IK 9 ' 9 . The results are presented in Figure 13 



Remark. Due to the vast computational expense to use DNS we are unable to show 
DNS results. 



4.5 Long time example 



We finally show a problem of the same form as example 4.2.1 but we will solve it 
for T = 0{e~ 2 ). In 1 120 II it was shown that the effective equation in this long time 
regime is of the form, 

t -Au xx -l3e 2 u xxxx = 0, Yx{0<t<T}, 
= f, u t = 0, Y x{t = 0}. 



*xxx ' 



This is still on the same flux form as assumed in Q with F = Au x + [3e 2 u x 
Therefore, it turns out that we only need to make the HMM process a little bit 
more accurate for long time computations. The modifications needed are: 

• Initial data in micro solver needs to be of higher order. We use a third order 
polynomial to approximate the higher macro derivatives. 

• The integration kernel needs to be smoother to give more accurate F (error 
less than 0[e 2 )) in order to capture the correct dispersion relationship, i.e. 
[e/r]) q < e 2 . This implies also that: 
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3D HOM (XY plane) 3D HMM (XY plane) 




0.5 1 0.5 1 



Figure 13: Three dimensional solutions of the homogenized equation and by using 
the HMM technique. 
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Figure 14: ID longtime DNS simulation (thin line) compared to a finite difference 
solution of the effective equation ([67]) (circles) and a HMM solution (crosses). 



The micro box needs to be a little bigger, t, 17 ~ e 1 2i lq , where q is defined in 
d2Ji 



We present the numerical computations in Figure 14 



5 Conclusions 

We have developed and analyzed numerical methods for multi-scale wave equa- 
tions with oscillatory coefficients. The methods are based on the framework of the 
heterogeneous multi-scale method (HMM) and have substantially lower compu- 
tational complexity than standard discretization algorithms. Convergence proofs 
for finite time approximation are presented in the case of periodic coefficients in 
multiple dimensions. Numerical experiments in one, two and three spatial dimen- 
sions show the accuracy and efficiency of the new techniques. Finally we explored 
simulation over very long time intervals. The effective equation for very long time 
is different from the finite time homogenized equation. Dispersive effects enter, 
and the effective equation must be modified 1 1 20 II . It is interesting to note that our 
HMM approach with just minor modifications accurately captures these dispersive 
phenomena. 



A Numerical schemes 

We present a detailed description of the numerical schemes used in the macro and 
micro solvers. The schemes are designed for one, two, three dimensions and can 
be generalized to higher dimensions. All the schemes are second order accurate in 
both time and space. 
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A.1 ID equation 

The finite difference scheme on the macro level has the form 



Fl 



(68) 



V 1 m±l/2 



: F (X m ±l/2> -Pmim)' 



where P n (u n - U n ,1 and P" ,, = i f [/" -[/"). The micro level 

m-l/2 H V m m—1 J m+1/2 H V m +l m 7 

scheme defined analogously: 



(u n+l 



2u n -u n - l + k 2 y n , 

mm J m 3 

m ~ T (/m+1/2 _ /m-l/2) ' 



J m+1/2 
, J m-l/2 



(69) 



ft 

u n —u 



■ a„ 1 ■ 



m-l 



A.2 2D equation 

The two dimensional problem is discretized with a scheme with the following 
schemes: The finite difference scheme on the macro level 



(70) 




where P^ + 1 ^ is given by (see Figure 



2.11 



P" , 

m+-e. 



L7 m _ ci +[/„_ 
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and the other P^ ± j ^ are components defined analogously. The micro level scheme 
is formulated as 

' ,.n+l _ 9 n _ n-1 , .2 n 



Y n = l( f m _ f m )+l(fW _ f (2) A 



m+ie. 



f (D 
J m-ie, 



(12) 



m—e 2 m+e 



,(n) 



2h 

(12) 



IZfi j 



ft 

,(21) 



a l 

-{<-< , ) + — r 1 ! 



iz" + iz" 

m+e 2 ^ u m-e 1 +e 2 



r(2) _ 
J m+\e 2 



a , ,,n _l_,,n 



a 



2fr 

(21) 



( 



+e. 



m— e, 1 m- 



iz" + lZ n 

m—e 2 m-e l -e 2 



(22) 



) 



r(2) 
\ 3 m-\e 2 



'» "1 /' ii m+e 1 ~ SrU m+e 1 -e 2 



2h 



2 



a 



h 

(22) 



7 ft V m m " e 2^ 

(72) 

When approximating /^_ le we take the average of u n m±e ^ and uJJ, +ei±f , 2 to approx- 
imate u(x m+ i ei±e2 ,t"). Then we use those two averages to approximate the y 
derivate of u at u(x m _i c ). The scheme is second order in both space and time. 

A.3 3D equation 

The macro scheme for the three dimensional problem is of the form 
I U" = 2U" - U"- 1 + K 2 Y" 

m mm m J 



yn 



H 



p(l,n) _ p {l,n) 



+ - 



H 



p(.2,n) _ p(2,n) 



+ - 



H 



7 (3,n) 



— F 



(3,n) 



where P" , is defined as, 

m+|e 3 



(73) 



P n , 

m+-e 3 



2H V 



l+«3 ^ 



1 f Um+e 2 + ^m+e 2 +e 3 ^m-e 2 +^m-e 2 +e 3 A 

2H V 2 2 J 

J} [ U m+e 3 ~ U m) 



(74) 
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and the other P n 1 defined analogously. The micro level scheme is a second order 



accurate scheme defined analogous with the 2D scheme 

f,.n+l _ n„n _ n-l , r.2 n 

V " = lf f m _ f m V-ff (2) -f (2) ^1 + if f C3) -f (3} 1 

(11) (12) 
(1) m +2 e i f n _ n\ m + 2 e i ( U m+e l+ e 2 ^ U m+e 2 m+e 1 -e 2 m-e 2 "\ 

^m+ie, ft ^"m+ej ""J + 2 fl I 2 2 J 



2/1 

a ( 13 ) 

m+ip, ^ U L„ _!_„ "I" U „J_„ U 1_„ _„ "I" U 



m +2 ei /* ni+e 1 +e 3 ni+e 3 m-\-e 1 —e 3 m+e 3 



/ m+e 1 +e 3 "1+63 m+e 1 — e 3 ni+e 3 \ 

v 2 2 J 



a j a '1 ;i n +7," 7," +71" 

(!) _ m ~2 e i f n „ "\ m ~2 e i f m+e, ^"m- ei +e 2 U m-e 2 ^ U m- ei -e 2 \ 
J --^i" h V"- U m-eJ+ 2h { 2 2 J 

(13) 

. m 2 e i f m+e 3 m-ei+eg _ m-e 3 m-e 1 -e 3 A 

2/i V 2 2 J 

(21) (22) 
a , ,,n 4-7/" 7/ n 4-7/" a 1 

(2) m+-e 2 A " m+£l+f , 2 -t-" m+ei m-e^e; tn-ej \ rn+-e 2 f n _ 

/ '"+^2 2/1 V 2 2 J h V U '"+ e 2 M 

(23) 

a j ,,n 1 ,,n ..n i ,.n 

m+-e 2 /'"• m+e2+e3 ' m+e 3 m+e 2 -e 3 T m-e 3 "\ 

+ 2fl V 2 2 J 

(21) (22) 

f (2) _ m--e 2 r u m+ei ^ "m+ ei - e2 ^ "m- ei -e 2 A m--e 2 ^ n r \ 

J m-\e 2 - 2h V 2 2 J ft V 1 " 1 " m " e J 

(23) 

a 1 4-7," 71" 4-7," 

f u m+e 3 ^ u m-e 2 +e 3 _ U m-e 2 ^ U m-e 2 -e 3 \ 

2h V 2 2 J 

(31) 

a ,u" _ +u" u" _ +u" 



r(3) _ m + 2 e 3 f m+e 1 +e 3 T "m+ej "m-e 1 +e 3 T "m-ej A 

^+ie 3 " 2d 1 2 2 J 

(32) (33) 
a , „n J-ji" 77" 4- 77 n 0- 1 

"■+-e 3 ^ " m+e2+e3 -t-" m+e2 " m - e2+e3 I - " m -e 2 -S m+-e 3 /■ n ^ n >. 

2hV 2 2 J h v m + e 3 U ™J 

(31) 

(3) m--e 3 /- " m+ei ^ "m+ ei -e 3 m-e 1 m-e 1 -e 3 \ 

■>m-ie 3 2h V 2 2 J 

(32) (33) 
a ! ,.n 4. j.n „n 1 n a i 

2d v 2 2 Jfrl m m " e 3j 



(75) 
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